#*************************###
"""                       ###
    本文件用于读取nc、h5文件
    将源数据读取、进行匹配等操作
"""                       ###
#*************************###
import os,sys
import datetime
import shutil

import numpy as np

import netCDF4
import h5py
from glob2 import glob
from file_utils import *

def read_iquam_universal(in_file, y_, m_, d_=-1, h_=-1):
    """
    TODO: 从.nc文件中 根据不同的条件 读取浮标数据
    :param in_file: 输入文件的路径
    :param y_: 所需数据的年份
    :param m_: 所需数据的月份
    :param d_: 所需数据的日期; if d_ != -1 该条件被激活
    :param h_: 所需数据的小时; if h_ != -1 该条件被激活
    :return obs_info: 一个包含浮标数据的字典
    """

    with netCDF4.Dataset(in_file) as nf:

        obs_info = {}
        obs_time = nf.variables["time"][:].data
        obs_year = nf.variables["year"][:].data
        obs_month = nf.variables["month"][:].data
        obs_day = nf.variables["day"][:].data
        obs_hour = nf.variables["hour"][:].data
        obs_ql = nf.variables["quality_level"][:].data
        obs_depth = nf.variables["depth"][:].data
        obs_type = nf.variables["platform_type"][:].data

        if d_ == -1:  # d_为-1代表”日“未激活，”日“、”时“都不参与筛选
            condition_idx = np.where(np.logical_and.reduce(
                (obs_year == y_, obs_month == m_, obs_depth < 10)))
        elif h_ == -1:  # h_为-1代表”时“未激活，“时”不参与筛选
            condition_idx = np.where(np.logical_and.reduce(
                (obs_year == y_, obs_month == m_, obs_day == d_, obs_depth < 10)))
        else:
            condition_idx = np.where(np.logical_and.reduce(
                (obs_year == y_, obs_month == m_, obs_day == d_, obs_hour == h_, obs_depth < 10)))

        ##********** time **********
        obs_info['year'] = nf.variables["year"][:][condition_idx]
        obs_info['month'] = nf.variables["month"][:][condition_idx]
        obs_info['day'] = nf.variables["day"][:][condition_idx].data
        obs_info['hour'] = nf.variables["hour"][:][condition_idx].data
        obs_info['mins'] = nf.variables["minute"][:][condition_idx].data
        obs_info['sec'] = nf.variables["second"][:][condition_idx].data
        obs_info['time'] = nf.variables["time"][:][condition_idx].data

        ##********** 大气压强 **********
        obs_info['air_pres'] = nf.variables["air_pressure"][:][condition_idx].data

        ##********** 大气温度 **********
        obs_info['air_temp'] = nf.variables["air_temperature"][:][condition_idx].data

        ##********** 云层含水量 **********
        obs_info['cc'] = nf.variables["cloud_coverage"][:][condition_idx].data

        ##********** 海温标记必须为0 **********
        obs_info['flag'] = nf.variables["sst_flags"][:][condition_idx].data

        ##********** 浮标观察质量 **********
        obs_info['ql'] = nf.variables["quality_level"][:][condition_idx].data

        ##********** 经纬度、深度 **********
        obs_info['dep'] = nf.variables["depth"][:][condition_idx].data
        obs_info['lat'] = nf.variables["lat"][:][condition_idx].data
        obs_info['lon'] = nf.variables["lon"][:][condition_idx].data

        ##********** 洋面温度 ************
        obs_info['sst'] = nf.variables["sst"][:][condition_idx].data

    return obs_info


def read_iquam_time_range(i_file, range_1, range_2):
    """
    按时间范围读取浮标数据
    :param i_file: 输入文件路径
    :param range_1: 时间范围界限1
    :param range_2: 时间范围界限2
    :return:
    obs_flag：浮标标记
    obs_sst：浮标海温
    obs_lat：纬度
    obs_lon, 经度
    obs_ql[condition_idx]：浮标质量
    obs_depth[condition_idx]：浮标深度
    obs_type[condition_idx]：浮标类型
    obs_time[condition_idx]：时间
    """
    with netCDF4.Dataset(i_file) as nf:
        obs_time = nf.variables["time"][:]
        obs_ql = nf.variables["quality_level"][:]
        obs_depth = nf.variables["depth"][:]
        obs_type = nf.variables["platform_type"][:]
        condition_idx = np.where(np.logical_and.reduce(
            (os.logical_and(obs_time >= range_1, obs_time <= range_2), obs_depth < 10, obs_ql > 4)))
        obs_mins = nf.variables["minute"][:][condition_idx]
        obs_flag = nf.variables["sst_flags"][:][condition_idx]
        obs_sst = nf.variables["sst"][:][condition_idx]
        obs_lat = nf.variables["lat"][:][condition_idx]
        obs_lon = nf.variables["lon"][:][condition_idx]

    return obs_flag, obs_sst, obs_lat, obs_lon, obs_ql[condition_idx], obs_depth[condition_idx], obs_type[
        condition_idx], obs_time[condition_idx]


def read_hy2b_l2b_sst(in_file, land_ocean_flag_=0):
    """
    TODO: 从h5文件当中读取 hy2b卫星l2b产品数据
    :param in_file:输入文件路径
    :param sat_flag_: 数据类型（为0就是海洋）
    :return sat_info: 一个包含卫星数据的字典
    """
    sat_file = h5py.File(in_file, 'r')
    sat_info = {}

    sat_info['land_ocean_flag'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Land_Ocean_Flag"].value
    shape_of_mat = sat_info['land_ocean_flag'].shape
    sat_sst_t = sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Res0_SST_Retrieve_Quality"].value

    condition_index = np.where(
        np.logical_and(np.logical_and(sat_info['land_ocean_flag'] == land_ocean_flag_, sat_sst_t != -9999.0),
                       np.logical_and(sat_sst_t != -8888.0, sat_sst_t != -7777.0)))

    sat_info['land_ocean_flag'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Land_Ocean_Flag"].value[condition_index]
    sat_info['abnormanl_flag'] = \
        np.zeros_like(sat_info['land_ocean_flag'])

    # ***************** sst ****************#
    ########################################
    sat_info['sst'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Res0_SST"].value[condition_index]
    sat_info['sst_ql'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Res0_SST_Retrieve_Quality"].value[condition_index]

    # **************** space ***************#
    ########################################
    sat_info['lat'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Lat_Of_Product"].value[condition_index]
    sat_info['lat'] = np.divide(sat_info['lat'], 1e6)
    sat_info['lon'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Long_Of_Product"].value[condition_index]
    sat_info['lon'] = np.divide(sat_info['lon'], 1e6)
    # sat_info['offNadir'] = sat_file["production_information/OffNadir"]

    # **************** time ****************#
    ########################################
    time_trans = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Scan_Time_Trans"].value

    sat_info['time'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Scan_Time"].value

    sat_info['orbit_direction'] = sat_file["foundation_information"].attrs['OrbitDirection'][0]

    sat_info['day_night_flag'] = -1  # 默认异常
    if sat_info['orbit_direction'] == b'ASCENDING':
        sat_info['day_night_flag'] = 0  # 夜晚
    elif sat_info['orbit_direction'] == b'DESCENDING':
        sat_info['day_night_flag'] = 1  # 白天

    sat_info['year'] = (time_trans[:, 0].reshape(-1, 1) * np.ones(shape_of_mat))[condition_index]
    sat_info['month'] = (time_trans[:, 1].reshape(-1, 1) * np.ones(shape_of_mat))[condition_index]
    sat_info['day'] = (time_trans[:, 2].reshape(-1, 1) * np.ones(shape_of_mat))[condition_index]
    sat_info['hour'] = (time_trans[:, 3].reshape(-1, 1) * np.ones(shape_of_mat))[condition_index]
    sat_info['min'] = (time_trans[:, 4].reshape(-1, 1) * np.ones(shape_of_mat))[condition_index]
    sat_info['sec'] = (time_trans[:, 5].reshape(-1, 1) * np.ones(shape_of_mat))[condition_index]

    # **************** other ***************#
    #########################################
    # 降雨标识 "0-NoRain,1-Rain"
    sat_info['rain_flag'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Rain_Flag"].value[condition_index]
    # 海冰标识 "0-Ocean,1-Ice"
    sat_info['ice_flag'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Ice_Flag"].value[condition_index]

    # 降雨率及其观测质量
    sat_info['ap'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Res0_AP"].value[condition_index]
    sat_info['ap_ql'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Res0_AP_Retrieve_Quality"].value[condition_index]

    # 云液水含量及其观测质量
    sat_info['cl'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Res0_CL"].value[condition_index]
    sat_info['cl_ql'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Res0_CL_Retrieve_Quality"].value[condition_index]

    # 海冰密集度及其观测质量
    sat_info['ic'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Res0_IC"].value[condition_index]
    sat_info['ic_ql'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Res0_IC_Retrieve_Quality"].value[condition_index]

    # 洋面风速、高风速及其观测质量
    sat_info['ssw'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Res0_SSW"].value[condition_index]
    sat_info['hssw'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Res0_hSSW"].value[condition_index]
    sat_info['ssw_ql'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Res0_SSW_Retrieve_Quality"].value[condition_index]

    # 大气水汽含量及其观测质量
    sat_info['wv'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Res0_WV"].value[condition_index]
    sat_info['wv_ql'] = \
        sat_file["data_fields/Res0_Retrieve_Swath_Fast_Product/Res0_WV_Retrieve_Quality"].value[condition_index]

    sat_file.close()

    return sat_info

def read_hy2b_l2d_sst(in_file):
    """
    TODO:解析h5py的hy2b_l2d文件
    :param in_file: 输入文件路径
    :return sat_info: 包含卫星数据的字典，每个字段用numpy数组存储
    """sat_info = {}
    sat_file = h5py.File(in_file, 'r')
    sat_sst_t = sat_file["data_fields/SST_SG"].value
    sat_ap_t = sat_file["data_fields/AP_SG"].value
    sat_cl_t = sat_file["data_fields/CL_SG"].value
    sat_ic_t = sat_file["data_fields/IC_SG"].value
    sat_wv_t = sat_file["data_fields/WV_SG"].value

    '''********************  选出海温数据正常的点  ********************'''
    mask = sat_sst_t != -9999.0
    mask = np.logical_and(mask, sat_sst_t != -8888.0)
    mask = np.logical_and(mask, sat_sst_t != -5555.0)
    '''
    mask = np.logical_and(mask, sat_ap_t != -9999.0)
    mask = np.logical_and(mask, sat_ap_t != -8888.0)
    mask = np.logical_and(mask, sat_ap_t != -5555.0)

    #mask = np.logical_and(mask, sat_cl_t != -9999.0)
    #mask = np.logical_and(mask, sat_cl_t != -8888.0)
    mask = np.logical_and(mask, sat_cl_t != -5555.0)

    #mask = np.logical_and(mask, sat_ic_t != -9999.0)
    #mask = np.logical_and(mask, sat_ic_t != -8888.0)
    mask = np.logical_and(mask, sat_ic_t != -5555.0)

    mask = np.logical_and(mask, sat_wv_t != -9999.0)
    mask = np.logical_and(mask, sat_wv_t != -8888.0)
    mask = np.logical_and(mask, sat_wv_t != -5555.0)
    '''

    condition_index = np.where(mask)

    # ***************** sst ****************#
    ########################################
    sat_info['sst'] = \
        sat_file["data_fields/SST_SG"].value[condition_index]
    sat_info['sst'] = sat_info['sst']/100 + 273.15
    sat_info['sst_ql'] = \
        sat_file["data_fields/SST_SG_Retrieve_Quality"].value[condition_index]

    # **************** space ***************#
    ########################################
    lat_t = sat_file["data_fields/Lat_SG_0.25"].value
    lon_t = sat_file["data_fields/Long_SG_0.25"].value
    lat_size = lat_t.shape[0]
    lon_size = lon_t.shape[1]
    lat_t = lat_t*np.ones((lat_size,lon_size))
    lon_t = lon_t*np.ones((lat_size,lon_size))
    sat_info['lat'] = \
        lat_t[condition_index]
    sat_info['lat'] = np.divide(sat_info['lat'], 1e3)

    sat_info['lon'] = \
        lon_t[condition_index]
    sat_info['lon'] = np.divide(sat_info['lon'], 1e3)
    # sat_info['offNadir'] = sat_file["production_information/OffNadir"]

    # **************** time ****************#
    ########################################

    sat_info['day_time'] = sat_file["data_fields/Scan_Time"].value[condition_index]
    sat_info['orbit_direction'] = sat_file["foundation_information"].attrs['OrbitDirection'][0]

    sat_info['day_night_flag'] = -1  # 默认异常
    if sat_info['orbit_direction'] == b'ASCENDING':
        sat_info['day_night_flag'] = 0  # 夜晚
    elif sat_info['orbit_direction'] == b'DESCENDING':
        sat_info['day_night_flag'] = 1  # 白天

    # **************** other ***************#
    #########################################

    # 降雨率及其观测质量
    sat_info['ap'] = \
        sat_file["data_fields/AP_SG"].value[condition_index]
    sat_info['ap_ql'] = \
        sat_file["data_fields/AP_SG_Retrieve_Quality"].value[condition_index]

    # 云液水含量及其观测质量
    sat_info['cl'] = \
        sat_file["data_fields/CL_SG"].value[condition_index]
    sat_info['cl_ql'] = \
        sat_file["data_fields/CL_SG_Retrieve_Quality"].value[condition_index]

    # 海冰密集度及其观测质量
    sat_info['ic'] = \
        sat_file["data_fields/IC_SG"].value[condition_index]
    sat_info['ic_ql'] = \
        sat_file["data_fields/IC_SG_Retrieve_Quality"].value[condition_index]

    # 洋面风速、高风速及其观测质量
    sat_info['ssw'] = \
        sat_file["data_fields/SSW_SG"].value[condition_index]
    sat_info['ssw_ql'] = \
        sat_file["data_fields/SSW_SG_Retrieve_Quality"].value[condition_index]

    # 大气水汽含量及其观测质量
    sat_info['wv'] = \
        sat_file["data_fields/WV_SG"].value[condition_index]
    sat_info['wv_ql'] = \
        sat_file["data_fields/WV_SG_Retrieve_Quality"].value[condition_index]

    sat_file.close()

    return sat_info

def time_and_space_match_l2b(
        out_txt, #输出txt的路径
        sat_dir, #卫星文件（h5）的路径
        obs_dir, #浮标文件（nc）的路径
        year_,   #观测的年份
        month_, #观测的月份
        day_=-1,#观测的日期
        hour_=-1,#观测的小时
        time_thr=10, #时间阈值
        space_thr=0.1, #空间阈值
):
    """
    将卫星数据与浮标数据进行匹配并输出到txt
    :param out_txt: 输出txt文档的路径
    :param sat_dir: 输入卫星数据文件路径
    :param obs_dir: 输入浮标数据文件路径
    :param year_: 需要匹配的年份
    :param month_: 需要匹配的月份
    :param day_: 需要匹配的日期
    :param hour_: 需要匹配的小时
    :param time_thr: 时间分辨率窗口阈值
    :param space_thr: 空间分辨率窗口阈值
    :return: None
    """
    ##两个阈值，
    # 时间匹配窗口阈值time_thr,
    # 空间匹配窗口阈值space_thr，
    # 这里设定为卫星观测的尺度，一个或两个像元，
    # 这里设为一个像元0.1度
    #time_thr = 10
    #space_thr = 0.01

    ###这里设定匹配的输出结果，也可以不落地
    #out_txt=r"F:\PIE\sst_blend\osisaf_avhrr_sst_iquam_20210519_20210526.txt"
    compare_info=open(out_txt,'a+') #追加写入txt

    ##从nc中读取浮标数据
    #obs_info中存储一个字典，
    #浮标中的数据以ndarray的格式存储在其各个字段中
    obs_info = read_iquam_universal(obs_dir,year_,month_,day_,hour_)

    ##从h5中读取卫星数据
    # sat_info中存储一个字典，
    # 浮标中的数据以ndarray的格式存储在其各个字段中
    sat_info = read_hy2b_l2d_sst(sat_dir)
    sat_sample_num = len(sat_info['sst'])

    count = 0
    for i in range(sat_sample_num):
        ##设置卫星观测的时间窗口，根据阈值
        time_str = '%d%02d%02d%02d%02d%02d' % (sat_info['year'][i],
                                               sat_info['month'][i],
                                               sat_info['day'][i],
                                               sat_info['hour'][i],
                                               sat_info['min'][i],
                                               sat_info['sec'][i])

        ref_time_range_1 = datetime.datetime.strptime(time_str,"%Y%m%d%H%M%S") - datetime.timedelta(minutes=time_thr)
        ref_time_range_2 = datetime.datetime.strptime(time_str,"%Y%m%d%H%M%S") + datetime.timedelta(minutes=time_thr)

        ##为方便跟iquam的时间进行比较，将时间转换为相同参考时间
        delta_time_range_1 = (ref_time_range_1 - datetime.datetime.strptime("19810101000000", "%Y%m%d%H%M%S")).total_seconds()
        delta_time_range_2 = (ref_time_range_2 - datetime.datetime.strptime("19810101000000", "%Y%m%d%H%M%S")).total_seconds()

        ##先进行时间匹配，（空间上iquam是覆盖全球先进行空间匹配无意义）
        obs_time_condition = np.logical_and(
            obs_info['time'] >= delta_time_range_1, obs_info['time'] <= delta_time_range_2)

        ##设置卫星经纬度窗口
        lat_s = sat_info['lat'][i]-space_thr
        lat_n = sat_info['lat'][i]+space_thr
        lon_w = sat_info['lon'][i]-space_thr
        lon_e = sat_info['lon'][i]+space_thr

        ##设置一个大概的空间范围来筛选可匹配的iquam，sst_s,sst_e,sst_n,sst_w是根据当前数据集中的经纬度四至坐标的最值来确定
        obs_space_condition = np.logical_and(np.logical_and(obs_info['lat'] >= lat_s, obs_info['lat'] <= lat_n),
                                             np.logical_and(obs_info['lon'] >= lon_w, obs_info['lon'] <= lon_e))


        ##将时空两个条件进行筛选过滤
        obs_idx = np.where(np.logical_and(obs_time_condition, obs_space_condition))

        compare_nums = len(obs_idx[0])
        if compare_nums > 0:
            compare_lat = obs_info['lat'][obs_idx[0]]
            compare_lon = obs_info['lon'][obs_idx[0]]
            count+=1
        else:
            continue

        # 存储各点与卫的星经纬度差
        difflat = compare_lat - sat_info['lat'][i]
        difflon = compare_lon - sat_info['lon'][i]

        # 存储各点与卫星的欧氏距离
        diff = (difflat ** 2) + (difflon ** 2)
        diff = diff ** 0.5

        #print(count)
        sample_min = np.argmin(diff) #最小距离的点在筛选出来点列表中的下标
        idx_min = (obs_idx[0][sample_min])#[0] #最小距离的点在整体样本中的下标
        print(year_,month_,day_,'len:',len(obs_idx),sample_min,idx_min,obs_idx[0])
        if str(obs_info['sst'][idx_min]) == "--":
            continue

        ##将匹配数据的各字段写入txt
        compare_line = [

            ##***************************************##
            #***************** 卫星信息 ****************#
            str(sat_info['sst'][i]/100+273.15), str(sat_info['sst_ql'][i]),  # 海温及其观测质量 0,1
            str(sat_info['land_ocean_flag'][i]), #陆地海洋标记（必为0） 2
            str(sat_info['abnormanl_flag'][i]), #异常标记（必为0） 3
            str(sat_info['year'][i]),str(sat_info['month'][i]),str(sat_info['day'][i]),#年月日 4,5,6
            str(sat_info['hour'][i]), str(sat_info['min'][i]), str(sat_info['sec'][i]),#时分秒 7,8,9
            str(sat_info['lat'][i]), str(sat_info['lon'][i]),#经纬度 10,11

            str(sat_info['rain_flag'][i]), str(sat_info['ice_flag'][i]),  # 降雨标记、海冰标记 12,13
            str(sat_info['ap'][i]), str(sat_info['ap_ql'][i]),  # 降雨率及其观测质量 14,15
            str(sat_info['cl'][i]), str(sat_info['cl_ql'][i]),  # 海冰密度及其观测质量 16,17
            str(sat_info['ssw'][i]), str(sat_info['hssw'][i]), str(sat_info['ssw_ql'][i]),  # 风速、高风速，风速观测质量 18,19,20
            str(sat_info['wv'][i]), str(sat_info['wv_ql'][i]),  # 大气水汽含量及其观测质量 21,22

            str(sat_info['orbit_direction']),str(sat_info['day_night_flag']), #23,24
            #str(sat_info['offNadir']),

            ##***************************************##
            # ***************** 浮标信息 ***************#

            # 25,26,27
            str(obs_info['year'][idx_min]), str(obs_info['month'][idx_min]), str(obs_info['day'][idx_min]),  # 年月日 4,5,6
            # 28,29,30
            str(obs_info['hour'][idx_min]), str(obs_info['mins'][idx_min]), str(obs_info['sec'][idx_min]),  # 时分秒 7,8,9
            # 31
            str(obs_info['time'][idx_min]),

            str(obs_info['dep'][idx_min]),  # 深度 32
            str(obs_info['lat'][idx_min]),  # 经度 33
            str(obs_info['lon'][idx_min]),  # 维度 34

            #
            str(obs_info['air_pres'][idx_min]),#气压 35
            str(obs_info['air_temp'][idx_min]),#气温 36
            str(obs_info['cc'][idx_min]),      #云层厚度 37
            str(obs_info['flag'][idx_min]),    #必须为0 否则不是海面温度 38
            str(obs_info['ql'][idx_min]),      #观测质量1-5 39

            str(obs_info['sst'][idx_min])]   #洋面温度 40

        compare_line_info = ",".join(compare_line) + "\n"
        compare_info.write(compare_line_info)
        #compare_info.close()

def time_and_space_match_l2d(
        out_txt, #输出txt的路径
        sat_dir, #卫星文件（h5）的路径
        obs_dir, #浮标文件（nc）的路径
        year_,   #观测的年份
        month_, #观测的月份
        day_=-1,#观测的日期
        hour_=-1,#观测的小时
        time_thr=10, #时间阈值
        space_thr=0.1, #空间阈值
        file_num = -1,
        cur_file = -1,
        data_count = -1
):
    """
    将卫星数据与浮标数据进行匹配并输出到txt
    :param out_txt: 输出txt文档的路径
    :param sat_dir: 输入卫星数据文件路径
    :param obs_dir: 输入浮标数据文件路径
    :param year_: 需要匹配的年份
    :param month_: 需要匹配的月份
    :param day_: 需要匹配的日期
    :param hour_: 需要匹配的小时
    :param time_thr: 时间分辨率窗口阈值
    :param space_thr: 空间分辨率窗口阈值
    :param file_num: 对应日期有多少h5文件需要匹配
    :param cur_file: 当前已经提取到多少日期
    :param data_count: 在txt中已写入多少条数据
    :return: None
    """
    ##两个阈值，
    # 时间匹配窗口阈值time_thr,
    # 这里选择30分钟或60分钟，代表时间窗口为正负30或60分钟，
    # 空间匹配窗口阈值space_thr，
    # 这里设定为卫星观测的尺度，一个或两个像元，
    # 这里设为一个像元0.25度
    #time_thr = 30
    #space_thr = 0.25

    ###这里设定匹配的输出结果，也可以不落地
    #out_txt=r"F:\PIE\sst_blend\osisaf_avhrr_sst_iquam_20210519_20210526.txt"
    compare_info=open(out_txt,'a+') #追加写入txt

    ##从nc中读取浮标数据
    #obs_info中存储一个字典，
    #浮标中的数据以ndarray的格式存储在其各个字段中
    obs_info = read_iquam_universal(obs_dir,year_,month_,day_,hour_)

    ##从h5中读取卫星数据
    # sat_info中存储一个字典，
    # 浮标中的数据以ndarray的格式存储在其各个字段中
    sat_info = read_hy2b_l2d_sst(sat_dir)
    sat_sample_num = len(sat_info['sst'])

    count = data_count
    for i in range(sat_sample_num):

        ##设置卫星观测的时间窗口，根据阈值
        if sat_info['day_time'][i]>24:
            continue
        sat_day_index = sat_dir.rfind('T')-2
        sat_day = int(sat_dir[sat_day_index:sat_day_index+2])
        sat_hour = int(sat_info['day_time'][i])
        sat_mins = int((sat_info['day_time'][i]-sat_hour)*60)
        sat_sec  = int((sat_info['day_time'][i]-sat_hour-sat_mins/60)*3600)

        time_str = '%04d%02d%02d%02d%02d%02d' % (year_,
                                               month_,
                                               sat_day,
                                               sat_hour,
                                               sat_mins,
                                               sat_sec)
        #t("sec",sat_sec,"day_time",sat_info['day_time'][i])
        #print("time_str",time_str)
        ref_time_range_1 = datetime.datetime.strptime(time_str,"%Y%m%d%H%M%S") - datetime.timedelta(minutes=time_thr)
        ref_time_range_2 = datetime.datetime.strptime(time_str,"%Y%m%d%H%M%S") + datetime.timedelta(minutes=time_thr)

        ##为方便跟iquam的时间进行比较，将时间转换为相同参考时间
        delta_time_range_1 = (ref_time_range_1 - datetime.datetime.strptime("19810101000000", "%Y%m%d%H%M%S")).total_seconds()
        delta_time_range_2 = (ref_time_range_2 - datetime.datetime.strptime("19810101000000", "%Y%m%d%H%M%S")).total_seconds()

        ##先进行时间匹配，（空间上iquam是覆盖全球先进行空间匹配无意义）
        obs_time_condition = np.logical_and(
            obs_info['time'] >= delta_time_range_1, obs_info['time'] <= delta_time_range_2)

        ##设置卫星经纬度窗口
        lat_s = sat_info['lat'][i]-space_thr
        lat_n = sat_info['lat'][i]+space_thr
        lon_w = sat_info['lon'][i]-space_thr
        lon_e = sat_info['lon'][i]+space_thr

        ##设置一个大概的空间范围来筛选可匹配的iquam，sst_s,sst_e,sst_n,sst_w是根据当前数据集中的经纬度四至坐标的最值来确定
        obs_space_condition = np.logical_and(np.logical_and(obs_info['lat'] >= lat_s, obs_info['lat'] <= lat_n),
                                             np.logical_and(obs_info['lon'] >= lon_w, obs_info['lon'] <= lon_e))

        ##将时空两个条件进行筛选过滤
        obs_idx = np.where(np.logical_and(obs_time_condition, obs_space_condition))

        ##输出交互
        print(year_, month_, sat_day,
              'cur_file: [%03d/%03d]'% (cur_file,file_num),
              '#'+'='*int((cur_file/file_num)*15)+'>'+'-'*(15-int(cur_file/file_num*15))+'#',
              "  cur_data: [%05d/%05d]" % (i, sat_sample_num),
              '#'+'='*int((i/sat_sample_num)*15)+'>'+ '-'*(15-int(i/sat_sample_num*15))+'#',
              "  matched:%06d"%count,)
        compare_nums = len(obs_idx[0])
        if compare_nums > 0:
            compare_lat = obs_info['lat'][obs_idx[0]]
            compare_lon = obs_info['lon'][obs_idx[0]]
            count+=1
        else:
            continue

        # 存储各点与卫的星经纬度差
        difflat = compare_lat - sat_info['lat'][i]
        difflon = compare_lon - sat_info['lon'][i]

        # 存储各点与卫星的欧氏距离
        diff = (difflat ** 2) + (difflon ** 2)
        diff = diff ** 0.5

        #print(count)
        sample_min = np.argmin(diff) #最小距离的点在筛选出来点列表中的下标
        idx_min = (obs_idx[0][sample_min])#最小距离的点在整体样本中的下标
        print(year_,month_,day_,'len:',len(obs_idx),sample_min,idx_min,obs_idx[0])
        if str(obs_info['sst'][idx_min]) == "--":
            continue

        ##将匹配数据的各字段写入txt
        compare_line = [

            ##***************************************##
            #***************** 卫星信息 ****************#
            #0,1
            str(sat_info['sst'][i]), str(sat_info['sst_ql'][i]),  # 海温及其观测质量
            #2,3,4
            str(year_),str(month_),str(sat_day),#年月日
            #5,6,7
            str(sat_hour), str(sat_mins), str(sat_sec),#时分秒
            #8
            str(sat_info['day_time'][i]),
            #9,10
            str(sat_info['lat'][i]), str(sat_info['lon'][i]),#经纬度


            #11,12
            str(sat_info['ap'][i]), str(sat_info['ap_ql'][i]),  # 降雨率及其观测质量
            #13,14
            str(sat_info['cl'][i]), str(sat_info['cl_ql'][i]),  # 海冰密度及其观测质量
            #15,16
            str(sat_info['ssw'][i]), str(sat_info['ssw_ql'][i]),  # 风速、高风速，风速观测质量
            #17,18
            str(sat_info['wv'][i]), str(sat_info['wv_ql'][i]),  # 大气水汽含量及其观测质量

            #19,20
            str(sat_info['orbit_direction']),str(sat_info['day_night_flag']),

            ##***************************************##
            # ***************** 浮标信息 ***************#

            # 21,22,23
            str(obs_info['year'][idx_min]), str(obs_info['month'][idx_min]), str(obs_info['day'][idx_min]),  # 年月日 4,5,6
            # 24,25,26
            str(obs_info['hour'][idx_min]), str(obs_info['mins'][idx_min]), str(obs_info['sec'][idx_min]),  # 时分秒 7,8,9
            # 27
            str(obs_info['time'][idx_min]),

            str(obs_info['dep'][idx_min]),  # 深度 28
            str(obs_info['lat'][idx_min]),  # 经度 29
            str(obs_info['lon'][idx_min]),  # 维度 30

            #
            str(obs_info['air_pres'][idx_min]),#气压 31
            str(obs_info['air_temp'][idx_min]),#气温 32
            str(obs_info['cc'][idx_min]),      #云层厚度 33
            str(obs_info['flag'][idx_min]),    #必须为0 否则不是海面温度 34
            str(obs_info['ql'][idx_min]),      #观测质量1-5 35

            str(obs_info['sst'][idx_min])]   #洋面温度 36

        compare_line_info = ",".join(compare_line) + "\n"
        compare_info.write(compare_line_info)

    compare_info.close()

def select_data_from_day_range(
        out_txt,year_,month_,day1_,day2_,sat_dir,obs_dir,time_thr=10,space_thr=0.1):
    """
    选取特定时间范围段进行时空匹配
    :param out_txt: 输出路径
    :param year_: 年份
    :param month_: 月份
    :param day1_:  开始日期
    :param day2_:  结束日期
    :param sat_dir:  卫星数据文件
    :param obs_dir:  浮标数据文件
    :return: None
    """
    for day in range(day1_,day2_+1):
        sat_paths = []
        #获取符合时间范围的卫星、浮标数据文件路径
        sat_paths = glob(sat_dir+'\\H2B_OPER_SMR_L2D_SG_%04d%02d%02dT??????_%04d%02d%02dT??????_???_????_??.zip'%(
            year_,month_,day,year_,month_,day
        ))

        sat_paths+=glob(sat_dir+'\\H2B_OPER_SMR_L2D_SG_%04d%02d%02dT??????_%04d%02d%02dT??????_???_????_??.tar.gz'%(
            year_,month_,day,year_,month_,day
        ))

        obs_paths = []
        obs_paths = glob(obs_dir+'\\%04d%02d-STAR-L2i_GHRSST-SST-iQuam*.nc'%(
            year_,month_
        ))

        f_count = 0

        for obs_path in obs_paths :
            for sat_comp in sat_paths: #sat_comp为压缩文件所在位置
                print("cur_file:",f_count)
                f_count+=1
                data_count = len(open(out_txt, 'rU').readlines())
                sat_super_path = '' #解压缩到的文件夹名称
                sat_mid_path = '' #可能产生的中间文件路径
                sat_path = '' #解压完成后的数据文件路径

                if sat_comp.find(".zip")!=-1:
                    sat_super_path = un_zip(sat_comp)
                    sat_path = sat_super_path+"\\"+sat_comp[sat_comp.rfind("\\")+1:].replace(".zip",".h5")
                elif sat_comp.find(".tar.gz")!=-1:
                    sat_mid_path = un_gz(sat_comp)
                    sat_super_path = un_tar(sat_mid_path)
                    sat_path = sat_super_path+"\\"+sat_comp[sat_comp.rfind("\\")+1:].replace(".tar.gz",".h5")

                time_and_space_match_l2d(out_txt,
                                     sat_path,
                                     obs_path,
                                     year_,
                                     month_,
                                     day,
                                     cur_file=f_count,
                                     file_num=len(sat_paths),
                                     data_count=data_count,
                                     space_thr=space_thr,
                                     time_thr=time_thr)

                if sat_comp.find(".zip")!=-1:
                    shutil.rmtree(sat_super_path)
                elif sat_comp.find(".tar.gz")!=-1:
                    shutil.rmtree(sat_super_path)
                    os.remove(sat_mid_path)





